Valence bond solid phases in a cubic antiferromagnet 



O 
O 

(N 

O 

Q 

in 



■4— » 
•4-* 



o 
o 



> 

(N 
(N 

O 

-i— > 
c3 



c 

o 
o 



X 



K. S. D. Beach and Anders W. Sandvik 
Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 02215 

(Dated: December 5, 2006) 

We report on a valence bond projector Monte Carlo simulation of the cubic lattice quantum 
Heisenberg model with additional higher-order exchange interactions in each unit cell. The model 
supports two different valence bond solid ground states. In one of these states, the dimer pattern is 
a three-dimensional analogue of the columnar pattern familiar from two dimensions. In the other, 
the dimers are regularly arranged along the four main diagonals in 1/8 of the unit cells. The phases 
are separated from one other and from a Neel phase by strongly first order boundaries. Our results 
strengthen the case for exotic transitions in two dimensions, where no discontinuities have been 
detected at the Heisenberg Neel- VBS transition driven by four-spin plaquet interactions. 



In dimension greater than one, quantum antiferromag- 
nets with nonfrustrating, local interactions are generi- 
cally Neel ordered at zero temperature. It is possible, 
however, for competing interactions to stabilize other 
singlet ground states that have only short range mag- 
netic correlations. These include valence bond solids 
(VBS) which spontaneously break the translational 
symmetry of the lattice, and spin liquids 0, which are 
featureless states having no broken symmetries. 

Spin liquids came to prominence following Anderson's 
proposal that a resonating valence bond (RVB) descrip- 
tion of the doped Mott insulator might provide a route to 
superconductivity in the cuprates [2|, |3j ■ We have since 
learned that spin liquids are quite unusual states of mat- 
ter. Whereas Neel and VBS phases have spin-1 excita- 
tions, spin liquids have spin-^ (spinon) excitations Q and 
topological order [f| . The existence of liquid phases has 
been established in quantum dimer models, but for phys- 
ical spin models there is at most suggestive evidence [fj. 
Because of Monte Carlo sign problems, none of the pur- 
ported spin liquid models is amenable to exact simulation 
except on very small clusters. 

Given the limited progress that has been made, it 
is worthwhile to build a repertoire of sign-problem-free 
models that exhibit quantum disordered phases. In two 
dimensions (2D), a good candidate is the square-lattice 
Heisenberg model with a ring exchange term that per- 
mutes the spins around an elementary plaquet. It is 
known that ring exchange can destroy magnetic order 
and drive the system into a VBS phase, both for XY 
and fully SU(2)-invariant [H models. It has been argued 
that the transition between antiferromagnetism and VBS 
order occurs at a special quantum critical point with de- 
confined spinon excitations 0]. Such a transition may 
also take place between two VBS phases with different 
dimer orderings [Id, EH- if this scenario holds, then it 
may be possible — having first indentified a deconfincd 
critical point — to move along some new axis in phase 
space (corresponding to some additional interaction) and 
follow the line of deconfined quantum critical points in 
the hope that it eventually opens into an extended spin 
liquid phase. As it turns out, the transitions in the XY 



case are either weakly first order 13, 13] or possibly con- 
tinuous but not of the type suggested by Senthil et al. 
141 ]. In the SU(2)-symmetric case, however, recent work 
provides compelling evidence of deconfined quantum crit- 
icality 0. 

The situation in three dimensions (3D) is less well un- 
derstood. Bernier and coworkers have studied the quan- 
tum Heisenberg model on the cubic lattice with nearest- 
and next-nearest-neighbour interactions [lq . Using an 
Sp(iV) generalization of the spin algebra, they show that 
there are three stable spins liquid phases for large N. 
These phases have short range magnetic correlations at 
wavevectors (0,0, tt), (0,tt, tt), and (tt, tt, 7t). It is un- 
likely that these phases survive in the physical N = 1 
limit. In particular, it appears certain that the (0,0, tt) 
and (tt, 7T, tt) states undergo a confinement transition to a 
VBS state. Because of the nature of the precursor mag- 
netic correlations, the translational symmetry breaking 
occurs either in one lattice direction, leading to a 3D 
generalization of the columnar phase familiar from 2D, 
or it occurs in all three directions, leading to a "box" 
VBS. Montrunich and Senthil have studied the Mott in- 
sulating phases of strongly correlated bosons in 3D and 
have reached similar conclusions about the possible va- 
lence bond ordering 17 1 . 

In this Letter we discuss two different 3D generaliza- 
tions of ring exchange. In both cases, the interaction is 
constructed as a product of four spin-spin operators in 
each unit cell. In addition to the Neel phase, we find two 
VBS phases with ordering vectors (0,0, tt) and (tt, tt,tt). 
Unlike the 2D case, the transitions between these phases 
are all strongly first order. We find no evidence of spin 
liquid behaviour. 

We consider a system of S = ^ spins arranged in a 
cubic lattice of even linear size L, subject to periodic 
boundary conditions. The Hamiltonian 

- H = J^P + U^PPPP + V^PPPP (1) 

H □ 

is expressed in terms of spin singlet projection operators 
Pij = -7 — Si ■ Sj whose site indices are understood to 
range over nearest neighbours (— ) and along the main 
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FIG. 1: The Hamiltonian contains interactions that act (a) 
between nearest neighbours, (b) along the diagonals of each 
cube, and (c) along the edges of each cube face. All interac- 
tions are built from singlet projection operators between spins 
in opposite sublattices, indicated by open and filled circles. 



diagonals (M) and face edges (□) of each cubic unit cell. 
The set of allowed interactions is illustrated in Fig. Q] 

The connection to physical spin operators is made by 
expanding the products (j — Sj ■ Sj)(4 — Sfc • S;) • ■ ■ that 
appear in Eq. ([1]). The first term is a two-spin interac- 
tion, whereas the second and third terms are alternating- 
sign linear combinations of two-, four-, six- and eight- 
spin interactions. We also considered an interaction 
—H ~ K PP acting around the cube faces (i.e., square 
plaqucts symmetrized in the three orthogonal directions), 
but as in the XY case [l8|, it is not sufficient to disrupt 
the Neel order — at least not in the K > parameter 
range that can be simulated. 

The simulations were carried out using the valence 
bond projector Monte Carlo algorithm [19]. In this 
scheme, the ground state is obtained by repeated ap- 
plication of the Hamiltonian to a singlet trial state: 
\tp) = lim T1 _ i . oo (-fl")™|'0 trial ). We selected as trial state 
a factorizable RVB wavefunction with bond amplitudes 
of the form h(r) = r~ p [20]. The exponent p was deter- 
mined self-consistently by measuring the powerlaw tail of 
the bond distribution in the final projected state. 

Our main result is that the system has three phases, 
characterized by a nonzero sublattice magnetization and 
by two forms of long range dimer order. In the thermo- 
dynamic limit, the operators 
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acquire a nonzero expectation value in each of the re- 
spective phases. Here, a — 1,2,3 is a label denoting 
the vector components and S = ei + e2 + e3 is a vector 
spanning the main diagonal of each cubic cell. In prac- 
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FIG. 2: Point styles identify data computed in the extreme 
J (•), U (■), and V (□) limits. (Top left) The VBS phases 
have dimer correlation functions T>\ = ^(D 81 • D H ) (solid 

lines) and D 2 = (D D • D D ) (dotted lines) very near sat- 
uration. (Top right) In the extreme J limit, the sublat- 
tice magnetization M — 0.4222(6) is within 1% of the spin 
wave theory result M sw = 0.42165. 22]. (Bottom left) The 
energy E/J = —1.652334(5) is also close to the predicted 
value E sv ,/J — —1.6455. In the large U limit, the energy 
E/U = -0.13053(6) is only about 4% below -1/8, which is 
the expectation value of the Hamiltonian in the static valence 
bond pattern. (Bottom right) The Neel state has gapless exci- 
tations. The extreme U state has a large gap A/U = 0.474(1) 
to triplet excitations. The extreme V state is also fully spin 
gapped, but the value of the gap cannot be reliably measured 
using operator-string resampling [2lj . 



tice, these quantities are measured on finite lattices of 
increasing size and extrapolated to the thermodynamic 
limit according to M 2 = lim^o^M-M), etc. See Fig. [2] 
The necessary loop estimators for the two- and four-spin 
operators are derived in Ref. |21| . 

In the extreme J limit (i.e., U, V — > 0), the Hamilto- 
nian reduces to the nearest neighbour quantum Heisen- 
berg model. The ground state is Neel ordered (M ^ 0) 
and has a continuous degeneracy corresponding to the 
spatial direction of the sublattice magnetization. The 
two VBS phases are connected to the large U and large 
V limits. In the one case, the translational symmetry 
is completely broken and there is an eightfold degener- 
ate ground state of the form D H ~ (±1, ±1, ±1). In the 
other, there is a sixfold degenerate ground state of the 
form D D - (±1, 0, 0), (0, ±1, 0), (0, 0, ±1), which breaks 
translational symmetry in only one of the three lattice 
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FIG. 3: A greyscale map of the average number of bonds 
completely contained in each cubic unit cell; the values range 
from (white) to 4 (black). The measurements are made with 
reference to a cell at the origin (indicated by an arrow) whose 
eight sites are connected by exactly four valence bonds. In 
this example, a portion of the L = 10 lattice is shown as a 
set of slices stacked vertically. The columns correspond to the 
three extreme cases J = 1, U = V = (left) , U = 1, J — V — 
(centre), and V = l, J = U=0 (right). 



directions. Our simulations indicate that the valence 
bond order is very strong: the extrapolated values are 
D m = 0.98(1) x (1,1,1) in the extreme U limit and 
D n = 0.957(2) x (1,0,0) in the extreme V limit. Af- 
ter equilibration, the configurations are essentially locked 
into a particular dimer pattern. The tunnelling time be- 
tween energetically equivalent configurations grows ex- 
ponentially in the system size and exceeds the lifetime of 
our longest simulations even for L = 6. 

Figure [3] shows the dimer patterns as they occur in the 
wavefunction itself. These are imaged as follows. For 
each bond configuration, we identify a unit cell whose 
eight sites are connected by exactly eight valence bonds. 
This cell is translated to the origin and a map is made 
of the remaining cells by counting the number of valence 
bonds that have both endpoints contained in the same 
cell. In the Neel state, the dimer correlations do not 
persist beyond one or two lattice spacings. The large U 
phase forms a spacing-2 superlattice in which 1/8 of the 
cells have valence bonds along all four main diagonals 
[Fig. Ufa)]; the large V phase consists of planar slabs of 
colinear bonds so that half of the cells have four bonds 
aligned in the preferred direction [Fig. St b)]. There are 
weak fluctuations about these static configurations. In 
the large U case, these are predominantly relocations of 
a bond cluster to one of the nearest neighbour interstitial 
cells [Fig. Etc)]. For L < 12, the fluctuations interact in 




FIG. 4: (a) The U interaction favours a superlattice dimer or- 
dering in which every eighth cell contains four valence bonds 
along the main diagonal. The translational symmetry is bro- 
ken in all three lattice directions, (b) The V interaction leads 
to a dimer configuration that breaks the translational sym- 
metry in one direction only, (c) Deep in the VBS phases, 
fluctuations in the bond order are primarily due to intersti- 
tial defects. 

such a way that they preferentially align in one direction 
(as evidence by the faint pattern seen in the second and 
fourth slices of the centre row of Fig. [3]). As the lattice 
size increases, this effect diminishes and the rotational 
symmetry of the cubic lattice is restored. 

The M, D®, and D D phases are separated by first or- 
der phase boundaries. The upper panel of Fig. [5] shows 
the magnetization during the transition from Neel to 
VBS. There are two paths shown. In one, the simula- 
tion is first performed for J = 1, U = and a typical 
configuration from this run is used as the starting con- 
figuration for a run at incrementally larger U. A sec- 
ond path begins at U = 1, J = with the results fed 
into new runs at incrementally larger J. There is strong 
hysteretic behaviour beginning at L = 6 that becomes 
worse as the system size increases. Near the transition, 
the Monte Carlo sampling is increasingly dominated by 
rare tunnelling events between configurations with dimcr- 
ordered, short range valence bonds and those with res- 
onating, long range bonds. The magnetization decreases 
from 0.42 to ~ 0.39 before collapsing abruptly. 

The bottom panel of Fig. [5] traces the dimer correla- 
tions through the transition from VBS to VBS. Here, 
strong hysteresis sets in at L — 8. The simulation has 
difficulty moving directly between the two incompatible 
dimer orders so it accomplishes the bond reconfiguration 
in a two-step process via some intermediate, metastablc 
bond configuration. The particular pathway differs de- 
pending on the direction across the transition, but the 
general behaviour appears to be that translational sym- 
metry breaking is preceded by a process of cell substitu- 
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tion. This can be seen by looking at the mixed correla- 
tion function (D m -D ). As U/V decreases from 00, edge 
cubes are replaced by diagonal cubes while the transla- 
tional symmetry remains broken at (0, 0, n) and the di- 
agonal cubes remain uncorrelated between the occupied 
slabs. The symmetry breaking transition occurs when 
the diagonal cubes finally organize themselves within and 
between the slabs. On the other hand, as U/V increases 
from 0, the occupied diagonal cubes are replaced by edge 
cubes of arbitrary direction, forming a 3D version of the 
plaquet phase. The symmetry restoration from (tt, tt, 7t) 
to (0, 0, tt) comes when the edge cubes lock in a common 
direction. 

In 2D, recent studies of a Neel-VBS transition in a 
Heisenberg model including four-spin interactions have 
shown evidence of deconfined quantum criticality. There 
are no signs of discontinuities, the extracted critical ex- 
ponents have reasonable values, and an emergent U(l) 
symmetry is seen explicitly [15[. Our results presented 
here show that the 3D Neel-VBS transition is very dif- 
ferent, with strongly first order boundaries between the 
phases. We do not observe any spin liquid phase. 
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hundred Monte Carlo time steps. (Bottom panel) The dimer 
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as a function of V/U for fixed J — 0. Here, the system size 
is L = 8. The cube edge correlator T>2 is computed from 
both the V/U = (■) and V/U = 00 (□) limits. The mixed 
correlation function V12 (A) is shown evolving from right to 
left. Its discontinuity is a reliable marker of the change in 
dimer order. The inset shows schematically the three phases 
separated by lines of first order transitions. 
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